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Summary 

It  is  well  known  that  solutions  of  the  backwards  heat  equation  can  be 
approximated  by  solutions  of  a pseudoparabolic  equation  depending  on  a small 
parameter  e.  The  emphasis  In  this  paper  is  on  the  mathematical  problems 
which  arise  In  approximating  solutions  to  Initial-boundary  value  problems  for 
this  perturbed  equation.  The  approximation  procedure  we  propose  is  obtained 
through  the  development  of  a potential  theory  for  pseudoparabolic  equations, 
the  asymptotic  evaluation  of  certain  contour  integrals,  and  results  based  on 
a theorem  of  Levin  in  the  theory  of  entire  functions. 

I . Introduction. 

The  problem  of  constructing  solutions  to  the  heat  equation  backwards  in 

time  is  dne  of  the  classical  improperly  posed  problems  in  partial  differential 

3 

equations.  Mathematically  the  problem  can  be  formulated  (in  R ) in  the  fol- 

3 3 

lowing  manner:  For  D a bounded  domain  in  R find  u(x,t),  jc  e R , such 
that 

AjU  ■ u(  in  D x (0,to)  (1-la) 

u - 0 on  3D  x (0,tQ)  (1.1b) 

u(£,to)  - iKjt)  in  D (1.1c) 

where  t(jt)  is  a prescribed  function.  As  is  well  known  ((15],  [16])  in 
general  no  solution  exists  to  (1.1a)  - (1.1c)  and  if  it  does  the  solution 
does  not  depend  continuously  on  the  data  $(j()  in  any  reasonable  norm.  Une 
approach  to  "solving"  (1.1a)  - (1.1c)  is  the  method  of  quasi-reversibility  as 


♦Professor  of  Mathematics,  University  of  Delaware,  Newark,  Delaware,  U.S.A. 


Approved  for  public  relep 
distribution  unlimited. 


79 


initially  developed  by  Lattes  and  Lions  (111]),  and  it  is  a version  of  this 
approach  that  we  wish  to  discuss  in  this  paper.  The  version  we  have  in  mind 
is  to  replace  the  problem  Cl -la)  - tl.lc)  by  the  modified  problem 


+ A,u  «0  in  D x (0,t  ) 

3 o 

(1.2a) 

0 on  3D  x (o,t  ) 
o 

(1.2b) 

to)  - *(£)  in  D 

(1.2c) 

where  e is  a small  positive  parameter  ([7],  [9]).  Three  obvious  questions 
immediately  present  themselves: 

1)  Does  a unique  solution  to  (1.2a)  - (1.2c)  exist  for  every  e > 0 
and  does  it  depend  continuously  on  the  data  $(£)? 

% 2)  As  e -*•  0 does  the  solution  of  (1.2a)  - (1.2c)  approach  the  solution 

of  (1.1a)  - (1.1c)  (if  it  exists!)? 

3)  What  constructive  methods  are  available  for  approximating  the 
solution  of  (1.2a)  - (1.2c)  for  e > 0? 

In  most  discussions  of  the  method  of  quasi-reversibility  attention  is  usually 
focused  on  questions  1)  and  2),  and  the  answer  to  3)  is  normally  "use  finite 
difference  approximations  or  partial  eigenfunction  expansions".  However  for 
multi-dimensional  problems,  large  time  intervals,  and  small  values  of  the 
parameter  e,  such  methods  are  in  general  rather  impractical.  It  is  in  fact 
with  this  difficulty  in  mind  that  we  have  chosen  to  use  the  third  order  per- 
turbed equation  (1.2a)  Instead  of  the  fourth  order  equation 


u 


ut  + AjU  ■ 0 


(1.3) 


as  proposed  by  Lattes  and  Lions  ([11]).  As  will  be  seen,  such  a choice  will 
enable  us  to  use  potential  theoretic  methods  and  asymptotic  analysis  to 
examine  solutions  of  (1.2a)  - (1.2c)  for  small  e,  thus  providing  a more  practi- 
cal method  for  constructing  approximations  to  solutions  of  the  backwards  heat 
equation  in  higher  dimensions  and  for  large  time  intervals. 

Having  tailored  our  model  (1.2a)  - (1.2c)  to  give  a satisfactory  answer 
to  question  3)  does  not  of  course  excuse  us  from  dealing  with  the  first  two 
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questions!  However,  the  question  of  existence,  uniqueness,  and  continuous 

dependence  on  the  data  Is  well  known  (c.f.  [1])  and  the  second  question  has 

been  answered  by  Ewing  in  the  following  theorem  ([7]): 

Theorem  1:  Let  u(^,t)  be  a solution  of  (1.1a),  (1.1b)  such  that 

| |u(^c,tQ)  - <)>(J6)||  j < ||u(j£,0)||  2 < M>  where  are  positive  constants, 

L L 

and  let  v(^c,t)  be  the  solution  of  (1.2a)  - (1.2c)  for  £ » [log  M/{]  *.  Then 
for  every  t > 0, 

I lu-vl I o i C(t)e 
L 

where  C(t)  does  not  depend  on  e. 

An  obvious  drawback  of  Ewing's  result  is  that  the  error  bound  depends 
logarithmically  on  6 instead  of  HBlder  continuously  as  in  the  stabilized 
quasi-reversibility  method  of  Millar  ([15]).  A slmlliar  problem  also 
arises  if  one  uses  quasi-reversibility  methods  in  conjunction  with  (1.3), 
or,  iu  a different  direction,  Tikhonov's  regularization  method  for  solving 
the  backwards  heat  equation  (c.f.  [8],  [18]).  In  our  case  this  means  that 
in  order  to  achieve  accurate  results  we  must  choose  the  parameter  e to  be 
quite  small  and  assume  that  the  data  at  t “ tQ  is  measured  with  a high 
degree  of  accuracy . 

Although  our  analysis  in  what  follows  is  presented  for  the  heat  equation 

3 2 

in  R , analogous  results  can  also  be  obtained  in  R (see  Section  IV) . 

II.  Potential  Theory  for  the  Pseudo-Heat  Equation. 

Equation  (1.2a)  is  a particular  example  of  an  equation  of  pseudo- 
parabolic  or  Sobolev  type  and  is  usually  referred  to  as  the  pseudo-heat 
equation.  Equations  of  pseudoparabolic  type  appear  in  a variety  of  areas  of 
application  and  have  been  the  object  of  a considerable  amount  of  attention  in 
recent  years.  For  information  concerning  this  class  of  equations  as  well  as 
an  extensive  blblograph',  a suggest  consulting  the  recent  book  by  Carroll  and 
Shuwalter  ([1]).  In  this  section  we  shall  develop  a potential  theory  for  the 
pseudo-heat  equation  ([3],  [5])  with  the  aim  of  using  these  results  in 


Section  III  to  construct  approximations  to  the  solution  of  Cl. 2a)  - (1.2c). 
We  shall  develop  this  potential  theory  for  the  "forward"  initial-boundary 


value  problem 


E A,u  - u + A,u  ■ 0 In  Dx(0,t  ) 
i t t 3 o 


(2.1a) 

(2.1b) 


u = f(j£,t)  on  3D  x (0,tQ)  (2.1b) 

u(£,0)  - 0 in  D (2.1c) 

where  f(jc,t)  is  assumed  to  be  continuously  differentiable  and  reserve  until 
Section  III  the  modifications  and  applications  of  these  results  to  the  study 
of  (1.2a)  - (1.2c). 

We  begin  by  defining  the  fundamental  solution  for  (2.1a)  by  ([3],  [5]) 


t_T)  ‘ " 77F  !' 


2.  _ . 

expl-wR  + M <‘t~~Ti]  du 
1 - e (D 


I uj 6 

/c 

where  R - |jt  - j^|  for  £ c R^  and  the  path  of  integration  is  a circle  of 

radius  A traversed  counterclockwise  about  the  point  a m — . If  we  now  assume 

/e 

that  D is  a bounded,  simply  connected  domain  with  Lyapunov  boundary  3D  and 
let  p(£,x)  denote  a continuous  density  defined  on  3D  x [0,t  ],  we  can  define 
a pseudo-heat  potential  by 


“<*•'>  ■ 1 1 L 


3^57  r(R*  t-T)dsdr 


where  v denotes  the  unit  normal  on  3D  pointing  into  D.  In  order  for  (2.3)  to 
be  a solution  of  (2.1a)  - (2.1c)  it  is  necessary  to  choose  p(£,r)  such  that 
(2.1b)  is  satisfied.  To  this  end  we  differentiate  both  sides  of  (2.3)  with 
respect  to  t and  let  * 3D.  Then  from  the  discontinuity  properties  of 
metaharmonic  potentials  we  arrive  at  an  integral  equation  for  p(£,t)  of  the 


e3/Z  If  • (l  + l + Jfei  + 


t3 


where  T is  a Fredholm  integral  operator  over  3D,  is  a Volterra  Integral 
operator  over  [0,t],  and  is  a Volterra  integral  operator  over  3D  x [ 0, t ] . 

The  kernels  of  these  Integral  operators  depend  on  r(R,t)  and  its  derivatives 
(c .f . [3] ) . Due  to  the  fact  that  e > 0 we  can  deduce  the  following  Theorem 
([3]): 


Theorem  2:  Let  C(<*D  x [0,tQ])  be  the  Banach  space  of  continuous  functions 
defined  over  3D  x [0,tQ]  with  respect  to  the  maximum  norm.  Then 
("t  + + fei  + ^2^  ^ exists  as  a bounded  linear  operator  on  C(3D  x [Ojt^])  . 

Hence  any  solution  of  (2.1a)  - (2.1c)  can  be  represented  in  the  form  (2.3)  for 
some  continuous  density  p(£,x). 

In  order  to  actually  construct  p(Jjj,t)  from  (2.4)  it  is  necessary  to 
determine  accurate  approximations  to  r(R,t)  ahd  its  derivatives  for  small 
values  of  c and  R ^ 0,  t >_  0.  Clearly  it  suffices  to  consider  the  function 

K(R,t)  “ " | exp  [ -ojR  + M 

l-cw 


6 


(2.5) 


exp(-t/e)  1 , R t 

L <f  exp  [-y  — + r-]dy  . 

re  »i  J /e  e(l-y^) 

|u-l|  - « 

The  evaluation  of  K(R,t)  can  be  conveniently  divided  into  three  separate 
cases  ([5]): 

1)  R-0,  t“0(l) . In  this  case  we  expand  exp  ( ) in  its  Taylor 

e(l-y  ) 

series  in  t and  Integrate  termwise  to  arrive  at 


K(0,t)  - 


t exp  (-  / ) 


.A 


♦(1/2.  2;  l/e) 


(2.6) 


where  *(a,c;z)  denotes  the  confluent  hypergeometrlc  function.  From  the 
asymptotic  behaviour  of  ♦ (a,c;z)  for  large  z (c.f.[6])  we  have  from  (2.6) 


that 
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K (0 


, r N (1/,)n(3/,)  n N+I] 

,t)  - E 2",  2n  (f)  + o (f) 

L“-°  1 


(2.7) 


where  (a)n  ' IfST1  • 

2)  R>0.  fO(e).  We  again  expand  exp  ( t . ) in  powers  of  t and  inte- 

E(l-U  ) 

grate  termwise  to  obtain 


K(R,t) 


2 exp(-t/  ) 


“ 1 t 0+1  R n+1/ 2 R 

n-0  (n+1)!n!  (2^  Kn+1/2  <vT) 


(2,8) 


where  Kn(z)  denotes  the  modified  Bessel  function.  For  t - 0(e),  (2.8) 

can  be  readily  approximated  by  truncating  the  series  and  applying  the 
Clenahaw-Luke  method  of  backward  recursion  ([13],  Section  11.8).  An  esti- 
mate  on  the  convergence  ratio  of  (2.8)  can  be  obtained  by  using  the  In- 
equality 


I (-7:)  2 K 1,  (^li1/,  r(n+1/,) 

2/T  '2  /i  2 2 


(2.9) 


3)  R>0 , t-O(l) . By  deforming  the  circle  |w  - — | » ( in  (2.5)  onto  the 

/e 

contour  pictured  below  In  Figure  1 


Figure  1 

can  rewrite  K(R,t)  In  the  form 


K(R,t) 


ri/F 


exp[ * + z2]  g (z)dz 

/t  E 


(2.10) 
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» 

where 

g£(z)  - exp 

By  using  ad  hoc  methods,  it  Is  now  possible  to  obtain  a complete  asymptotic 
expansion  for  K(R,t)  as  e 0 in  the  form 
2 

K(R.O  " -fZ  **P(“  If)  [l+d^fj+djOf)^  ...  + dn(|)n  + OCCf)1^1)  1 (2.12) 

where  the  coefficients  dj  are  expressible  in  terms  of  Hermlte  polynomials 
([5]).  In  particular 


— “ 

ez4 


(2.11) 


H„ 


2*44  “8  2/t 


(-t:)  + 3 h6(-S=> 


4J  ” 2/t 

with  similar  expressions  holding  for  the  higher  order  coefficients 


(2.13) 


In 


(2.13)  Hn(z)  denotes  Hermlte 's  polynomial. 


The  expansions  in  all  of  the  above  cases  may  be  differentiated  termwise. 


Ill • Approximation  of  Solutions  to  the  "Backwards"  Pseudo-Heat  Equation. 

We  now  discuss  the  problem  of  approximating  solutions  of  initial- 
boundary  value  problems  for  the  "backwards"  pseudo-heat  equation  defined 
by  (1.2a)  - (1.2c).  One  approach  proceeds  as  follows.  By  replacing  t by 
to~t,  and  using  the  Fourier  transform  or  Borel  transform  to  construct  a solu- 
tion to  the  pure  initial  value  problem,  we  can  reduce  problem  (1.2a)  - (1.2c) 
to  an  initial-boundary  value  problem  of  the  form 

e 4jUt  - ut  - AjU  - 0 in  D x (0,tQ)  (3.1a) 

u ■ f(jc,t)  on  3D  x (0,tQ)  (3.1b) 

u(g,0)  « 0 in  D.  (3.1c) 

Note  that  (3.1a)  - (3.1c)  differs  from  (2.1a)  - (2.1c)  only  in  a sign  change  in 
the  differential  equation.  Although  this  has  no  effect  on  the  well-posedness 
of  the  problem,  as  we  shall  see  it  unfortunately  has  serious  implications  on 


the  problem  of  constructing  approximate  solutions ' for  small  values  of  the 
(positive)  parameter  e.  This  is,  of  course,  not  surprising  for  the  limiting 
case  as  e -*•  0 is  now  an  improperly  posed  problem,  i.e.  the  backwards  heat 
equation.  If  we  follow  the  analysis  of  the  previous  section,  it  is  seen  that 
we  can  represent  the  solution  of  (3.1a)  - (3.1c)  in  the  form 

.2 


u(^,t) 


_1 

2ir 


if  f 

T Jo  J 3D 


p(£,t) 


3v3x 


r(R,x-t)dsdx 


(3.2) 


where  p(£,x)  is  determined  as  the  solution  of  an  integral  equation  of 
Fredholm-Volterra  type  and  the  notation  is  the  same  as  in  (2.3).  We  observe 
that  the  only  difference  between  (3.2)  and  (3.3)  is  that  the  argument  of  the 
fundamental  solution  is  x-t  instead  of  t-x.  This  change,  however,  unfor- 
tunately complicates  the  evaluation  of  r(R,t)  since  we  are  now  interested 
in  the  asymptotic  behaviour  of  r(R,t)  for  t < 0 instead  of  t > 0.  Except 
for  the  expansion  (2.8)  the  analysis  of  Section  II  breaks  down  and  we  are 
forced  into  adopting  an  ad  hoc  approach  yielding  expansions  of  T(R,t)  in 
terms  of  a series  of  generalized  hypergeometric  functions  or  a Bessel- 
Laguerre  series  (c.f.[4J).  Such  expansions  are  of  limited  value  in  the  evalua- 
tion of  r(R,t)  for  t _ 0(1)  and  hence  we  are  led  to  look  for  other  approxi- 
mation procedures. 

To  present  such  an  alternate  approach  we  return  to  the  formulation 
(1.2a)  - (1.2c).  We  first  note  that  if  Xj  is  a sequence 


lim 

J-h.  2 > 0 

j 

then  in  a sufficiently  small  complex  neighbourhood  of  [0,tQ]  we  can  approximate 

any  analytic  function  defined  in  this  neighbourhood  by  a finite  linear  com- 

-X2  t “ 

bination  of  functions  taken  from  the  set  (e  j TM8  result  follows 

from  a theorem  of  Levin  in  the  theory  of  entire  functions  and  the  reader  is 
referred  to  [12],  p.  219,  for  details.  From  the  results  of  [2]  we  can  now 
conclude  that  any  continuous  function  defined  in  D can  be  approximated 
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in  D by  a linear  combination  of  the  functions 

r'1'2  V*/2  (V>  p»<co*  <J-‘> 

where  J (z)  denotes  Bessel's  function  and  Pm(z)  an  associated  Legendre 
n n 

function.  We  observe  that 


nmj 


(j^.t) 


_1.  X 

r 2 J .1.  (X.r)  Pm(co8  0)  exp  [im$ * — t]  (3.5) 

n+  / - j n - 


1 + eX 


j 


is  a solution  of  (1.2a).  Although  the  set  {u  . ) is  a complete  set  of 

nmj 

solutions  when  e-0((2)),  for  e>0  the  order  of  the  equation  is  Increased 
and  hence  except  for  exceptional  circumstances  we  would  expect  that  extra 

functions  would  have  to  be  Introduced  to  supplement  the  above  set.  (An 

2 

exceptional  case  occurs  for  example  when  D is  a sphere  and  {X^}  are  the 
eigenvalues  of  the  Laplaclan  in  D.)  In  particular  from  the  Runge  approxima- 
tion property  for  pseudoparabolic  equations  ([17])  and  the  results  of  Section 
II  (approximating  the  pseudo-heat  potential  defined  on  the  cylinder 

S x [0,to]  3 D x lO,tQ]  by  numerical  quadrature)  we  have  the  following 

2 

Theorem,  where  completeness  is  with  respect  to  the  L norm  over  D x [0,t  ] . 

o 

Theorem  3:  Let  {X.}  be  as  in  (3.3)  and  let  {£.}.  . be  a dense  set  of 
j j“l 

points  on  a sphere  S 73  D.  Then  the  functions 


wnj 


T"  tJff  r(Uj  - «!•  t-T)dT 


-S 


2 (Xjr)  Pn(cos  0)  exp  [imifi t] 


1+eX 


l4jl 


j 


form  a complete  set  of  solutions  to  (1.2a)  defined  in 


where  v 
D x [0,toJ . 

Now  let  (T  } , denote  the  complete  set  of  functions  defined  in 

n n- 1 

Theorem  3.  In  order  to  find  an  approximate  solution  to  (1.2a)  - (1.2c) 


we  minimize  the  functional 
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C.  -♦<«>llt2/m  + H * Vn<fre>l 

■1  L (D)  n-1 


i 2 
l/(3D 


[0,toD 


(3.6) 


N 

and  define  the  approximate  solution  to  be  ref  (x,t).  The  functions  ¥ 

n-1  " n * n 

can  be  evaluated  by  using  existing  tables  for  the  computation  of  Bessel  and 
Legendre  functions  in  combination  with  the  results  of  Section  II  to  approxi- 
32 

mate  r |£  - x|,  t-r) . Due  to  the  fact  that  (1.2a)  - (1.2c)  is  well 

posed,  the  minimization  procedure  indicated  above  is  stable.  Following 

Gaponenko  ([10])  we  require  that  the  minimization  of  (3.6)  allow  not  only  for 

the  variation  of  the  coefficients  c but  also  for  the  basis  functions  ¥ 

n n 

themselves,  i.e.  by  considering  £ j , j = 1,  ...,  N,  and  , j - 1,  ...,  N,  to 
be  variables  as  well.  Thus  we  are  led  to  the  problem  of  finding  the  minimum 
of  a function  of  the  6 N variables  Re  cn>  lm  cn»  * (£^,  5^), 

where  1 < ] < N.  A related  minimization  problem  is  considered  in  [10],  to- 
gether with  its  numerical  implementation  by  the  conjugate  gradient  method,  and 
a similar  approach  for  elliptic  boundary  value  problems  has  also  been  discussed 
by  Mathon  ([14]) . 


2 

IV.  The  Backwards  Heat  Equation  in  R . 

The  analysis  of  the  previous  sections  can  be  extended  to  treat  the 

2 

backwards  heat  equation  in  R , and  in  this  section  we  briefly  outline  the 
necessary  modifications.  Theorem  1 is  valid  in  any  number  of  dimensions.  A 
potential  theory  for  the  two  dimensional  pseudo-heat  equation  can  be  developed 
exactly  as  in  Section  II  where  the  fundamental  solution  is  now  defined  by 


(4.1) 


with  K (z)  denoting 
o 


a modified  Bessel  function  of  order  zero. 


From  [6]  we 
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